### Alizade, Dancygier, Ditlmann 
### "National Penalties Reversed"
### Replication Code 
### Figure 1
### For questions, contact jalizade@princeton.edu

# empty environment
rm(list = ls())


#setwd("")
setwd("C:/Users/Jey/Dropbox/WZB/NaturalizationExperiment/Submission/JOP/replication_JOP/data")


# load necessary packages
library(readstata13)
library(lmtest)
library(sandwich)
library(ggplot2)


# load data set
dat <- read.dta13("data_experimental.dta")


### regressions ###

# center-right
right_mod <- lm(e1_response ~ e1_treat_turkish + e1_treat_dual + e1_treat_benefits, dat[dat$e1_leftpol==0,])

# center-left
left_mod <- lm(e1_response ~ e1_treat_turkish + e1_treat_dual + e1_treat_benefits, dat[dat$e1_leftpol==1,])


# combine both models to a list
mods <- list(right_mod, left_mod)


# estimate robust standard errrors
robust_mods <- lapply(mods, function(x) coeftest(x, vcov = vcovHC(x, "HC1")))


### figure ###

# convert list with robust SEs to a data frame
df <- data.frame(do.call(rbind, robust_mods))


# remove rows with intercepts and keep columns with estimates and SEs (we do not need t-statistics and p-values)
df <- df[-c(1,5), 1:2]

# rename columns
names(df) <- c("coef", "se")

# treatment indicator
df$treat <- rep(c("Turkish", "Dual", "Benefits"), 2)

# partisanship indicator
df$pol <- rep(c("Center-Right", "Center-Left"), each = 3)
df$pol <- factor(df$pol, levels = c("Center-Left", "Center-Right"))

# create figure
fig1 <- ggplot(df, aes(x = treat, y = coef), group = pol)
fig1 <- fig1 + geom_hline(yintercept = 0, linetype = "dashed")
fig1 <- fig1 + geom_errorbar(aes(ymax = coef + se, ymin = coef - se,
                                 group = pol, color = factor(pol)), size = 1.2, 
                             position = position_dodge(width = 0.5), width = 0)
fig1 <- fig1 + geom_errorbar(aes(ymax = coef + 1.96*se, ymin = coef - 1.96*se,
                                 group = pol, color = factor(pol)), size = 0.4, 
                             position = position_dodge(width = 0.5), width = 0.2)
fig1 <- fig1 + geom_point(aes(x = treat, y = coef, group = pol, shape = factor(pol), fill = factor(pol), color = factor(pol)), 
                          position = position_dodge(width = 0.5), size = 3)
fig1 <- fig1 + xlab("") + ylab("Treatment Effect")
fig1 <- fig1 + coord_flip()
fig1 <- fig1 + scale_y_continuous(breaks= seq(-0.05, 0.2, 0.05))
fig1 <- fig1 + scale_shape_manual(values = c(24, 22), name = "Politician\nPartisanship")
fig1 <- fig1 + scale_fill_manual(values = c("red", "black"), name = "Politician\nPartisanship")
fig1 <- fig1 + scale_color_manual(values = c("red", "black"), name = "Politician\nPartisanship")
fig1 <- fig1 + guides(shape = guide_legend(reverse = T), color = guide_legend(reverse = T), fill = guide_legend(reverse = T))
fig1 <- fig1 + theme_bw()
fig1 <- fig1 + theme(text = element_text(size=24), panel.grid = element_blank())


